#OS
Dat <-read.csv("TCGA_allclinical.csv",header = T,row.names = 1,check.names = F)
Dat <-Dat [,9:25]
colnames(Dat)[1:2]<-c("time","status")
Dat$gender <- as.numeric(factor(Dat$gender))
Dat$MetStatus <- ifelse(Dat$MetStatus=="Metastatic",1,0)
Dat$stage <- as.numeric(factor(Dat$stage))
Dat$histological_type <- as.numeric(factor(Dat$histological_type))
#dat$histological_type <- as.numeric(factor(dat$histological_type))
#dat$SCNA_Cluster <- as.numeric(factor(dat$SCNA_Cluster))
#$chromosome.3.status <- as.numeric(factor(dat$chromosome.3.status))

pbc<-Dat

#row.names(pbc) <-pbc[,1]
#pbc<-pbc[,-1]


pbc$died <- pbc$status==1

head(pbc)
library(rms)




dd<-datadist(pbc)
options(datadist="dd")
options(na.action="na.delete")
summary(pbc$MetStatus)

coxpbc<-cph(formula = Surv(time,died) ~  age +gender+ histological_type+MetStatus+DL_riskscore,data=pbc,x=T,y=T,surv = T,na.action=na.delete)  #,time.inc =2920

print(coxpbc)

surv<-Survival(coxpbc) 
surv3<-function(x) surv(1095,x)
surv4<-function(x) surv(1825,x)


x<-nomogram(coxpbc,fun = list(surv3,surv4),lp=T,
            funlabel = c('3-year survival Probability','5-year survival Probability'),
            maxscale = 100,fun.at=c(.001,.01,.05,seq(.1,.9,by=.1),.95,.99,.999))



pdf("nomogram.pdf",width = 12,height = 10)
plot(x, lplabel="Linear Predictor",
     xfrac=.35,varname.label=TRUE, varname.label.sep="=", ia.space=.2, 
     tck=NA, tcl=-0.20, lmgp=0.3,
     points.label='Points', total.points.label='Total Points',
     total.sep.page=FALSE, 
     cap.labels=FALSE,cex.var = 1.6,cex.axis = 1.05,lwd=5,
     label.every = 1,col.grid = gray(c(0.8, 0.95)))
dev.off()

library(nomogramEx)
nomogramEx(nomo=x,np=2,digit = 9)


#不喜欢默认的颜色，先设置几个颜色
mycol<-c("#A6CEE3","#1F78B4","#33adff","#2166AC")
names(mycol) = c("dencol","boxcocl","obscol","spkcol")
mycol<- as.list(mycol)

#install.packages("regplot","vioplot","sm","beanplot")
library(survival)
pbccox <- coxph(formula = Surv(time,died) ~  age +gender+ histological_type+MetStatus+DL_riskscore
                , data = pbc)




library(regplot)

pdf("nomogram_new.pdf")
regplot(pbccox,
        observation=pbc[2,],
        failtime = c(1095,1825), 
        dencol="#FFD700",boxcol="#FFD700",
        prfail = TRUE, #cox回归中需要TRUE
        showP = T, #是否展示统计学差异
        droplines = F,#观测2示例计分是否画线
         #用前面自己定义的颜色
         #根据统计学差异的显著性进行变量的排序
        interval="confidence") #展示观测的可信区间


library(survival)






f3<-cph(formula = Surv(time,status) ~age +gender+ histological_type+MetStatus+DL_riskscore,data=pbc,x=T,y=T,surv = T,na.action=na.delete,time.inc = 1095) 

#参数m=50表示每组50个样本进行重复计算
cal3<-calibrate(f3, cmethod="KM", method="boot",u=1095,m=15,B=1000) 

pdf("calibration_3y.pdf",width = 8,height = 8)
plot(cal3,
     lwd = 2,#error bar的粗细
     lty = 1,#error bar的类型，可以是0-6
     errbar.col = c("#2166AC"),#error bar的颜色
     xlim = c(0,1),ylim= c(0,1),
     xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)",
     cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=0.6) #字的大小
lines(cal3[,c('mean.predicted',"KM")], 
      type = 'b', #连线的类型，可以是"p","b","o"
      lwd = 2, #连线的粗细
      pch = 16, #点的形状，可以是0-20
      col = c("#2166AC")) #连线的颜色
mtext("")
box(lwd = 1) #边框粗细
abline(0,1,lty = 3, #对角线为虚线
       lwd = 2, #对角线的粗细
       col = c("#224444")#对角线的颜色
) 
dev.off()


f5<-cph(formula = Surv(time,status) ~age +gender+ histological_type+MetStatus+DL_riskscore,data=pbc,x=T,y=T,surv = T,na.action=na.delete,time.inc = 1425) 
cal5<-calibrate(f5, cmethod="KM", method="boot",u=1325,m=15,B=1000)


plot(cal5,
     lwd = 2,
     lty = 1,
     errbar.col = c("#B2182B"),
     xlim = c(0,1),ylim= c(0,1),
     xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)",
     col = c("#B2182B"),
     cex.lab=1.2,cex.axis=1, cex.main=1.2, cex.sub=0.6)
lines(cal5[,c('mean.predicted',"KM")],
      type= 'b',
      lwd = 2,
      col = c("#B2182B"),
      pch = 16)
mtext("")
box(lwd = 1)
abline(0,1,lty= 3,
       lwd = 2,
       col =c("#224444"))


pdf("calibration_compare.pdf",width = 4,height = 4)
plot(cal3,lwd = 2,lty = 0,errbar.col = c("#DAA520"),
     bty = "l", #只画左边和下边框
     xlim = c(0,1),ylim= c(0,1),
     xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)",
     col = c("#DAA520"),
     cex.lab=1.2,cex.axis=1, cex.main=1.2, cex.sub=0.6)
lines(cal3[,c('mean.predicted',"KM")],
      type = 'b', lwd = 1, col = c("#DAA520"), pch = 16)
mtext("")

plot(cal5,lwd = 2,lty = 0,errbar.col = c("#DC143C"),
     xlim = c(0,1),ylim= c(0,1),col = c("#DC143C"),add = T)
lines(cal5[,c('mean.predicted',"KM")],
      type = 'b', lwd = 1, col = c("#DC143C"), pch = 16)

abline(0,1, lwd = 2, lty = 3, col = c("#224444"))

legend("topleft", #图例的位置
       legend = c("3-year","5-year"), #图例文字
       col =c("#DAA520","#DC143C"), #图例线的颜色，与文字对应
       lwd = 2,#图例中线的粗细
       cex = 1.2,#图例字体大小
       bty = "n")#不显示图例边框
dev.off()


#

#OS
library(survival)

Srv = Surv(pbc$time, pbc$died)

#此处选择5年的时间节点，输入文件的time列的单位是天，5年是1825天。
#下面每两行计算1种cox模型的系数，后面将画图对比



coxmod1 = coxph(Srv ~ age +gender+ histological_type+MetStatus+DL_riskscore, data=pbc)
pbc$nomogram = c(1 - (summary(survfit(coxmod1,newdata=pbc), times=1825)$surv))


coxmod2 = coxph(Srv ~ stage, data=pbc)
pbc$stage = c(1 - (summary(survfit(coxmod2,newdata=pbc), times=1825)$surv))
coxmod3 = coxph(Srv ~ histological_type, data=pbc)
pbc$histological_type = c(1 - (summary(survfit(coxmod3,newdata=pbc), times=1825)$surv))


head(pbc)
source("stdca.R") #stdca.R文件位于当前文件夹

mod1<-stdca(data=pbc, outcome="status", ttoutcome="time", timepoint=1825, 
            predictors=c("nomogram","chromosome.3.status","SCNA_Cluster"), cmprsk=TRUE, smooth=TRUE, xstop=0.5,intervention="FALSE")




pdf("net_benefit_OS.pdf",width = 6,height = 6)
stdca(data=pbc, outcome="status", ttoutcome="time", timepoint=1825,
      predictors=c("nomogram","chromosome.3.status","SCNA_Cluster"), 
      cmprsk=TRUE, smooth=TRUE, 
      xstop=0.5,intervention="FALSE")
dev.off()

###
library(caret)
library(ggDCA)

library(rms)
m1 <- lrm(status ~  age +gender+ histological_type+MetStatus+DL_riskscore,pbc,maxit=1000)
m2 <- lrm(status ~  histological_type,pbc,maxit=1000)
m3 <- lrm(status ~  MetStatus,pbc,maxit=1000)


d_m1 <- dca(m1,m2,m3)


library(ggprism)

pdf("net_benefit_OS.pdf",width = 6,height = 4)
ggplot(d_m1,linetype =F,lwd = 1.0)+
   theme_classic()+  
   theme_prism(base_size =17)+

  scale_colour_prism(         
      palette = "prism_light",
      name = "Cylinders",
   label = c("nomogram", "histological_type","Metastasis","ALL", "None"))
dev.off()



##ROC

library(survivalROC)
roc1 <- survivalROC(Stime=pbc$time, 
                    status=pbc$status, 
                    marker = pbc$nomogram, 
                    predict.time =1*365, # 计算pred.time时刻的ROC，一般是五年生存
                    method="KM")

roc2 <- survivalROC(Stime=pbc$time, 
                    status=pbc$status, 
                    marker = pbc$nomogram, 
                    predict.time =3*365, # 计算pred.time时刻的ROC，一般是五年生存
                    method="KM")

roc3 <- survivalROC(Stime=pbc$time, 
                      status=pbc$status, 
                      marker = pbc$nomogram, 
                      predict.time =4*365, # 计算pred.time时刻的ROC，一般是五年生存
                      method="KM")
auc1<- roc1$AUC





#par(mfrow = c(1,2))

plot(roc1$FP, roc1$TP, type="l", xlim=c(0,1), ylim=c(0,1),col="#FFD700", # 画峰值时的ROC
     xlab="1-Specificity (FPR)", ylab="Sensitivity (TPR)",
     lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
lines(roc2$FP, roc2$TP, type="l", col="#FF69B4", xlim=c(0,1), ylim=c(0,1),# 画峰值时的ROC
      
      lwd = 2)
lines(roc3$FP, roc3$TP, type="l", col="#1E90FF", xlim=c(0,1), ylim=c(0,1),# 画峰值时的ROC
      
      lwd = 2)
text(0.6,0.6,paste0("1-year AUC = ",round(roc1$AUC,3)),cex=1,col="#FFD700")
text(0.6,0.4,paste0("3-year AUC = ",round(roc2$AUC,3)),cex=1,col="#FF69B4")
text(0.6,0.2,paste0("5-year AUC = ",round(roc3$AUC,3)),cex=1,col="#1E90FF")

lines(x=c(0,1),y=c(0,1),lwd=1.5,lty=2,col="grey40")





